#######################################################################################################
#######Minorities overlooked: Group-based power-sharing and the exclusion-amid-inclusion dilemma#######
#########################################Replication Script############################################
#######################################################################################################

library(plyr)
library(dplyr)
library(rms)
library(stargazer)
library(data.table)
library(matrixStats)
library(ggplot2)
library(ggpubr)

######1. Data import######
cpsd <- read.csv("****ENTER PATH TO DATA FILE****")


######2. Descriptives######

###2.1. Table 2: Table of maximum power-sharers###
table2 <- cpsd[c("country","gwid","year","group","other","mminority","ps1h_corp_mo","ps1h_corp_g_strength","ps1h_lib_g_strength","ps1h_lib_mo","status_no")]
table2$ps1h_mo <- pmax(table2$ps1h_corp_mo, table2$ps1h_lib_mo)
table2$ps1h_mo_period <- round(as.numeric(table2$ps1h_mo), digits = 1)
table2_lag <- table2
table2_lag$year <- table2_lag$year + 1
table2_lag$ps1h_mo_period_lag <- table2_lag$ps1h_mo_period
table2 <- left_join(table2,table2_lag[c("country","group","year","ps1h_mo_period_lag")],by=c("country","group","year"))
table2$ps1h_mo_period_lag <- ifelse(is.na(table2$ps1h_mo_period_lag),-999,table2$ps1h_mo_period_lag)
table2$ps1h_mo_period <- ifelse(is.na(table2$ps1h_mo_period),-999,table2$ps1h_mo_period)
table2$pschange <- ifelse(table2$ps1h_mo_period != table2$ps1h_mo_period_lag,1,0)
table2$pschange <- ifelse(is.na(table2$pschange),1,table2$pschange)
table2 <- table2 %>% group_by(country,group) %>% arrange(country,group,year) %>% mutate(period = cumsum(pschange))#dito
table2$ps1h_mo_period <- ifelse(table2$ps1h_mo_period == -999, NA, table2$ps1h_mo_period)
table2$status_no <- as.numeric(paste(table2$status_no))
#min, max for each period
table2 <- table2 %>% group_by(country,group,period) %>% arrange(country,group,period) %>% mutate(status_no = median(status_no))
table2 <- table2 %>% group_by(country,group,period) %>% arrange(country,group,period) %>% mutate(from = min(year))
table2 <- table2 %>% group_by(country,group,period) %>% arrange(country,group,period) %>% mutate(to = max(year))
table2 <- table2 %>% group_by(country,group,period) %>% arrange(country,group,period) %>% mutate(max_psh = max(ps1h_mo))
table2 <- table2 %>% group_by(country,group,period) %>% arrange(country,group,period) %>% mutate(min_self = min(ps1h_corp_g_strength))
table2$min_self <- round(as.numeric(table2$min_self),digits=2)
table2 <- table2 %>% group_by(country,group,period) %>% arrange(country,group,period) %>% mutate(max_self = max(ps1h_corp_g_strength))
table2$max_self <- round(as.numeric(table2$max_self),digits=2)
table2 <- table2 %>% group_by(country,group,period) %>% arrange(country,group,period) %>% mutate(min_lib = min(ps1h_lib_g_strength))
table2$min_lib <- round(as.numeric(table2$min_lib),digits=2)
table2 <- table2 %>% group_by(country,group,period) %>% arrange(country,group,period) %>% mutate(max_lib = max(ps1h_lib_g_strength))
table2$max_lib <- round(as.numeric(table2$max_lib),digits=2)
#putting them together
table2$time_period <- ifelse(table2$from != table2$to, paste(table2$from,"-",table2$to,sep=""),table2$from)
table2$ps_corp_self <- ifelse(table2$min_self != table2$max_self, paste(table2$min_self,"-",table2$max_self,sep=""),table2$max_self)
table2$ps_lib <- ifelse(table2$min_lib != table2$max_lib, paste(table2$min_lib,"-",table2$max_lib,sep=""),table2$max_lib)
table2 <- table2[order(-table2$ps1h_mo_period,table2$country,table2$period,table2$other, table2$group),]
table2 <- unique(table2[c("country","group","mminority","time_period","ps_corp_self","ps_lib","status_no")])
write.csv(table2,file="table2.csv")

###2.2. Figure 2: Degree of corporate power-sharing of micro-minorities in different institutional regimes###
descriptive <- cpsd
#Creating factor variable for the type of institutionalized corporate power-sharing faced by each micro-minority
descriptive$type_corp <- ifelse(descriptive$ps1h_corpg_0_mm == 0 & descriptive$ps1h_corpg_0_nm == 0, "No accommodation",NA)
descriptive$type_corp <- ifelse(descriptive$ps1h_corpg_0 == 0 & (descriptive$ps1h_corpg_0_mm == 1 | descriptive$ps1h_corpg_0_nm == 1), "Only accommodation\nof other groups",descriptive$type_corp)
descriptive$type_corp <- ifelse(descriptive$ps1h_corpg_0 == 1 & descriptive$ps1h_corpg_0_nm == 1, "Accommodated alongside\nmain minority segments",descriptive$type_corp)
descriptive$type_corp <- ifelse(descriptive$ps1h_corpg_0 == 1 & descriptive$ps1h_corpg_0_nm == 0, "Accommodated without\nmain minority segments",descriptive$type_corp)
descriptive$type_corp <- as.factor(descriptive$type_corp)
#subset to only micro-minorities and summing up number of micro-minorities in each type
descriptive <- subset(descriptive, mminority == 1)
freq_corp <- descriptive[c("mminority","type_corp","one")]
freq_corp <- freq_corp %>% group_by(mminority) %>% arrange(mminority) %>% mutate(count_total = sum(one))
freq_corp <- freq_corp %>% group_by(mminority, type_corp) %>% arrange(mminority, type_corp) %>% mutate(freq_type_corp = sum(one) / count_total)
freq_corp <- unique(freq_corp)
#Figure 2: Frequency plot
setwd("****ENTER PATH TO EXPORT GRAPHS TO****")
figure2 <- ggplot(data=freq_corp, aes(x=type_corp, y=freq_type_corp)) +
  geom_bar(stat="identity") +
  scale_x_discrete(limits=c("No accommodation","Accommodated without\nmain minority segments","Accommodated alongside\nmain minority segments","Only accommodation\nof other groups")) +
  theme(axis.text = element_text(size=8), text=element_text(family="Times")) +
  labs(x = "Type of corporate power-sharing (micro-minority)", y = "Frequency")
ggsave(file="figure2.jpg", figure2, width = 13.8, height = 5, units ="cm", dpi = 800)

###2.3. Figure 3: "Mean" status of micro-minorities in different institutional regimes###
#formula for statistical summary
min.mean.sd.max <- function(x) {
  r <- c(min(x), mean(x) - sd(x), mean(x), mean(x) + sd(x), max(x))
  names(r) <- c("ymin", "lower", "middle", "upper", "ymax")
  r
}
#subset only to the cases used in the quantitative analysis, i.e. countries that have at least one micro-minority whose political status is coded by EPR (there is no reliable data on the de-facto status of the added "other" groups)
descriptive2 <- subset(descriptive, mminority_exist_epr == 1)
descriptive2 <- data.frame(descriptive2)
figure3 <- ggplot(aes(y = status_no, x = factor(type_corp)), data = descriptive2) + stat_summary(fun.data = min.mean.sd.max, geom = "boxplot") + xlab("Type of corporate power-sharing (micro-minority)") + ylab("'Mean' power status") +
  scale_x_discrete(limits=c("No accommodation","Accommodated without\nmain minority segments","Accommodated alongside\nmain minority segments","Only accommodation\nof other groups")) + 
  theme(axis.text = element_text(size=8), text=element_text(family="Times"))
ggsave(file="figure3.jpg", figure3, width = 13.8, height = 5, units ="cm", dpi = 800)


######3. Statistical Analysis######
#(again) subset only to the cases used in the quantitative analysis, i.e. countries that have at least one micro-minority whose political status is coded by EPR (there is no reliable data on the de-facto status of the added "other" groups)
cpsd_epr <- subset(cpsd,mminority_exist_epr == 1)

###3.1. Preliminaries###
#formula to decrease model size of glm's
stripGlmLR = function(cm) {
  cm$data = c()
  attr(cm$terms,".Environment") = c()
  attr(cm$formula,".Environment") = c()
  cm
}
#control variables to be included
control_vars <- "mminority + tek_state + geo_conc + ongoing_grp + d10_victory_neg + minoritysum + democracy + loggdppc + logpop + year"# + log(gdppc) + log(pop) 
control_vars_inc <- c("mminority", "tek_state", "geo_conc", "ongoing_grp", "d10_victory_neg", "minoritysum", "democracy", "loggdppc", "logpop", "year")

###3.2. Models###

##a) Model 1: Country indices##
model1 <- cpsd_epr[c("gwid","country","group","gwgroupid","year","ps1h_corp_mo","ps1h_lib_mo","status_no","included","discriminated", control_vars_inc)]
model1$gwid <- as.factor(paste(model1$gwid))
model1 <- model1[complete.cases(model1),]
dd=datadist(model1)
options(datadist="dd")
m1.1 <- lrm(as.formula(paste("status_no ~ gwid + ps1h_corp_mo + ps1h_lib_mo + ",control_vars)),data=model1, x =T, y =T, tol=1e-9, maxit= 200)
m1.1 <- robcov(m1.1,model1$gwid)
cse_m1.1 <- data.frame(sqrt(diag(vcov(m1.1))))
m1.2 <- stripGlmLR(glm(as.formula(paste("included ~ gwid + ps1h_corp_mo + ps1h_lib_mo + ",control_vars)), family = binomial(link = "logit"), data=model1))
cse_m1.2 <- data.frame(cluster.se(m1.2,as.integer(model1$gwid))[, 2])
m1.3 <- stripGlmLR(glm(as.formula(paste("discriminated ~ gwid + ps1h_corp_mo + ps1h_lib_mo + ",control_vars)), family = binomial(link = "logit"), data=model1))
cse_m1.3 <- data.frame(cluster.se(m1.3,as.integer(model1$gwid))[, 2])

##b) Model 2: Country indices, interactions##
m2.1 <- lrm(as.formula(paste("status_no ~ gwid + ps1h_corp_mo * mminority + ps1h_lib_mo * mminority + ",control_vars)),data=model1, x =T, y =T, tol=1e-9, maxit= 200)
m2.1 <- robcov(m2.1,model1$gwid)
cse_m2.1 <- data.frame(sqrt(diag(vcov(m2.1))))
m2.2 <- stripGlmLR(glm(as.formula(paste("included ~ gwid + ps1h_corp_mo * mminority + ps1h_lib_mo * mminority + ",control_vars)), family = binomial(link = "logit"), data=model1))
cse_m2.2 <- data.frame(cluster.se(m2.2,as.integer(model1$gwid))[, 2])
m2.3 <- stripGlmLR(glm(as.formula(paste("discriminated ~ gwid + ps1h_corp_mo * mminority + ps1h_lib_mo * mminority + ",control_vars)), family = binomial(link = "logit"), data=model1))
cse_m2.3 <- data.frame(cluster.se(m2.3,as.integer(model1$gwid))[, 2])

##c) Model 3: Group corporate index##
model3 <- cpsd_epr[c("gwid","country","group","gwgroupid","ps1h_corp_g_strength","ps1h_lib_mo","status_no","included","discriminated", control_vars_inc)]
model3$gwid <- as.factor(paste(model3$gwid))
model3 <- model3[complete.cases(model3),]
dd=datadist(model3)
options(datadist="dd")
m3.1 <- lrm(as.formula(paste("status_no ~ gwid + ps1h_corp_g_strength + ps1h_lib_mo + ",control_vars)),data=model3, x =T, y =T, tol=1e-9, maxit= 200)
m3.1 <- robcov(m3.1,model3$gwid)
cse_m3.1 <- data.frame(sqrt(diag(vcov(m3.1))))
m3.2 <- stripGlmLR(glm(as.formula(paste("included ~ gwid + ps1h_corp_g_strength + ps1h_lib_mo + ",control_vars)), family = binomial(link = "logit"), data=model3))
cse_m3.2 <- data.frame(cluster.se(m3.2,as.integer(model3$gwid))[, 2])
m3.3 <- stripGlmLR(glm(as.formula(paste("discriminated ~ gwid + ps1h_corp_g_strength + ps1h_lib_mo + ",control_vars)), family = binomial(link = "logit"), data=model3))
cse_m3.3 <- data.frame(cluster.se(m3.3,as.integer(model3$gwid))[, 2])

##d) Model 4: Group corporate index, interactions##
m4.1 <- lrm(as.formula(paste("status_no ~  gwid + ps1h_corp_g_strength * mminority + ps1h_lib_mo * mminority + ",control_vars)),data=model3, x =T, y =T, tol=1e-9, maxit= 200)
m4.1 <- robcov(m4.1,model3$gwid)
cse_m4.1 <- data.frame(sqrt(diag(vcov(m4.1))))
m4.2 <- stripGlmLR(glm(as.formula(paste("included ~ gwid + ps1h_corp_g_strength * mminority + ps1h_lib_mo * mminority + ",control_vars)), family = binomial(link = "logit"), data=model3))
cse_m4.2 <- data.frame(cluster.se(m4.2,as.integer(model3$gwid))[, 2])
m4.3 <- stripGlmLR(glm(as.formula(paste("discriminated ~ gwid + ps1h_corp_g_strength * mminority + ps1h_lib_mo * mminority + ",control_vars)), family = binomial(link = "logit"), data=model3))
cse_m4.3 <- data.frame(cluster.se(m4.3,as.integer(model3$gwid))[, 2])

##e) Model 5:  Group corporate index, interactions, strength of corporate index of other minorities##
model4 <- cpsd_epr[c("gwid","country","group","gwgroupid","ps1h_corp_g_strength","ps1h_lib_g_strength","ps1h_lib_mo","ps1h_corp_mo_others","ps1h_lib_mo_others","status_no","included","discriminated", control_vars_inc)]
model4$gwid <- as.factor(paste(model4$gwid))
model4 <- model4[complete.cases(model4),]
dd=datadist(model4)
options(datadist="dd")
m5.1 <- lrm(as.formula(paste("status_no ~  gwid + ps1h_corp_g_strength * mminority + ps1h_lib_mo + ps1h_corp_mo_others + ",control_vars)),data=model4, x =T, y =T, tol=1e-9, maxit= 200)
m5.1 <- robcov(m5.1,model4$gwid)
cse_m5.1 <- data.frame(sqrt(diag(vcov(m5.1))))
m5.2 <- stripGlmLR(glm(as.formula(paste("included ~ gwid + ps1h_corp_g_strength * mminority + ps1h_lib_mo + ps1h_corp_mo_others + ",control_vars)), family = binomial(link = "logit"), data=model4))
cse_m5.2 <- data.frame(cluster.se(m5.2,as.integer(model4$gwid))[, 2])
m5.3 <- stripGlmLR(glm(as.formula(paste("discriminated ~ gwid + ps1h_corp_g_strength * mminority + ps1h_lib_mo + ps1h_corp_mo_others + ",control_vars)), family = binomial(link = "logit"), data=model4))
cse_m5.3 <- data.frame(cluster.se(m5.3,as.integer(model4$gwid))[, 2])

##f) Model 6: Group corporate index, interactions, strength of corporate index of other minorities (interaction)##
m6.1 <- lrm(as.formula(paste("status_no ~  gwid + ps1h_corp_g_strength * mminority + ps1h_lib_mo + ps1h_corp_mo_others * mminority + ",control_vars)),data=model4, x =T, y =T, tol=1e-9, maxit= 200)
m6.1 <- robcov(m6.1,model4$gwid)
cse_m6.1 <- data.frame(sqrt(diag(vcov(m6.1))))
m6.2 <- stripGlmLR(glm(as.formula(paste("included ~ gwid + ps1h_corp_g_strength * mminority + ps1h_lib_mo + ps1h_corp_mo_others* mminority + ",control_vars)), family = binomial(link = "logit"), data=model4))
cse_m6.2 <- data.frame(cluster.se(m6.2,as.integer(model4$gwid))[, 2])
m6.3 <- stripGlmLR(glm(as.formula(paste("discriminated ~ gwid + ps1h_corp_g_strength * mminority + ps1h_lib_mo + ps1h_corp_mo_others* mminority + ",control_vars)), family = binomial(link = "logit"), data=model4))
cse_m6.3 <- data.frame(cluster.se(m6.3,as.integer(model4$gwid))[, 2])

##g) Model 7 (appendix 4): Group corporate index, interaction, group liberal index, strength of corporate index of other minorities, strength of liberal index of other minorities##
m7.1 <- lrm(as.formula(paste("status_no ~  gwid + ps1h_corp_g_strength * mminority + ps1h_lib_g_strength + ps1h_corp_mo_others + ps1h_lib_mo_others + ",control_vars)),data=model4, x =T, y =T, tol=1e-9, maxit= 200)
m7.1 <- robcov(m7.1,model4$gwid)
cse_m7.1 <- data.frame(sqrt(diag(vcov(m7.1))))
m7.2 <- stripGlmLR(glm(as.formula(paste("included ~ gwid + ps1h_corp_g_strength * mminority + ps1h_lib_g_strength + ps1h_corp_mo_others + ps1h_lib_mo_others + ",control_vars)), family = binomial(link = "logit"), data=model4))
cse_m7.2 <- data.frame(cluster.se(m7.2,as.integer(model4$gwid))[, 2])
m7.3 <- stripGlmLR(glm(as.formula(paste("discriminated ~ gwid + ps1h_corp_g_strength * mminority + ps1h_lib_g_strength + ps1h_corp_mo_others + ps1h_lib_mo_others + ",control_vars)), family = binomial(link = "logit"), data=model4))
cse_m7.3 <- data.frame(cluster.se(m7.3,as.integer(model4$gwid))[, 2])

###3.3. Export of models for the paper###
#main models
stargazer(type="html", m1.1, m2.1, m3.1, m4.1, m5.1, m6.1, se=c(cse_m1.1, cse_m2.1, cse_m3.1, cse_m4.1,cse_m5.1,cse_m6.1), title="Table 3: The effects of de-jure power-sharing on minorities' de-facto political status", column.labels = c("Power Status", "Power Status", "Power Status", "Power Status", "Power Status", "Power Status"),omit=c("gwid"),  dep.var.labels.include=FALSE, add.lines=list(c("Country-Fixed Effects","Yes","Yes","Yes","Yes","Yes")), style = "ajps", notes.append = TRUE, notes.align = "l", notes = "Country-clustered errors in parentheses")
#main models: included (appendix 4)
stargazer(type="html", m1.2, m2.2, m3.2, m4.2, m5.2, m6.2, se=c(cse_m1.2, cse_m2.2, cse_m3.2, cse_m4.2,cse_m5.2,cse_m6.2), title="The effects of de-jure power-sharing on minorities' de-facto government inclusion", column.labels = c("Included", "Included", "Included", "Included", "Included", "Included"),omit=c("gwid"),  dep.var.labels.include=FALSE, add.lines=list(c("Country-Fixed Effects","Yes","Yes","Yes","Yes","Yes")), style = "ajps", notes.append = TRUE, notes.align = "l", notes = "Country-clustered errors in parentheses")
#main models: discriminated (appendix 4)
stargazer(type="html", m1.3, m2.3, m3.3, m4.3, m5.3, m6.3, se=c(cse_m1.3, cse_m2.3, cse_m3.3, cse_m4.3,cse_m5.3,cse_m6.3), title="The effects of de-jure power-sharing on minorities' de-facto discrimination", column.labels = c("Discriminated", "Discriminated", "Discriminated", "Discriminated", "Discriminated", "Discriminated"),omit=c("gwid"),  dep.var.labels.include=FALSE, add.lines=list(c("Country-Fixed Effects","Yes","Yes","Yes","Yes","Yes")), style = "ajps", notes.append = TRUE, notes.align = "l", notes = "Country-clustered errors in parentheses")
#further models: splitting of average liberal PS index (appendix 4)
stargazer(type="html", m7.1, m7.2, m7.3, se=c(cse_m7.1, cse_m7.2, cse_m7.3), title="The effects of de-jure power-sharing on minorities [split-up liberal power-sharing index]", column.labels = c("Power Status", "Included", "Discriminated"),omit=c("gwid"),  dep.var.labels.include=FALSE, add.lines=list(c("Country-Fixed Effects","Yes","Yes","Yes")), style = "ajps", notes.append = TRUE, notes.align = "l", notes = "Country-clustered errors in parentheses")

###3.4. odds ratio###
exp(coef(m5.1))#odds ratio; ps1h_corp_g: 271.3467; ps1h_corp_g * mm: 0.06322925; liberal: 5.594386; corp_others: 0.1274203; 

###3.5. Graphical test of parallel odds assumption (model 5.1)###
#cf https://stats.idre.ucla.edu/r/dae/ordinal-logistic-regression/
library(Hmisc)
sf <- function(y) {
  c('Y>=0' = qlogis(mean(y >= 0)),
    'Y>=1' = qlogis(mean(y >= 1)),
    'Y>=2' = qlogis(mean(y >= 2)),
    'Y>=3' = qlogis(mean(y >= 3)))
}
(s <- with(model4, summary(as.numeric(paste(status_no)) ~ gwid + mminority * ps1h_corp_g_strength + ps1h_lib_mo + ps1h_corp_mo_others + tek_state + geo_conc + ongoing_grp + d10_victory_neg + minoritysum + democracy + loggdppc + logpop + year, fun=sf)))#all, including gwid
(s_small <- with(model4, summary(as.numeric(paste(status_no)) ~ mminority * ps1h_corp_g_strength + ps1h_lib_mo + ps1h_corp_mo_others + tek_state + geo_conc + ongoing_grp + d10_victory_neg + minoritysum + democracy + loggdppc + logpop + year, fun=sf)))#everything except FE's
png(file="parallelodds.png",width=4000,height=14000,units="px",res=300)
op <- par(mar=c(6,6, 4, 2) + 0.1)
p <- plot(s, which=1:4, pch=1:3, xlab='logit', main=' ', xlim=c(-5.472995,2.730029))
p
dev.off()
png(file="parallelodds_small.png",width=4000,height=4000,units="px",res=300)
op <- par(mar=c(6,6, 4, 2) + 0.1)
p <- plot(s_small, which=1:4, pch=1:3, xlab='logit', main=' ', xlim=c(-5.472995,2.730029))
p
dev.off()

###3.6. Predicted probabilities###
##obtaining the input independent variable values for the predictions##
prediction_values <- unique(subset(cpsd_epr[c("country","gwid","group","year", "geo_conc", "d10_victory_neg", "tek_state", "minoritysum","ongoing_grp","loggdppc","logpop","democracy")],year==2013))
prediction_values <- prediction_values %>% group_by(gwid,year) %>% arrange(gwid,year) %>% mutate(ongoing_grp = median(ongoing_grp,na.rm=T))#PS institutions targeting each group within a given country year on average
prediction_values <- prediction_values %>% group_by(gwid,year) %>% arrange(gwid,year) %>% mutate(geo_conc = median(geo_conc,na.rm=T))#PS institutions targeting each group within a given country year on average
prediction_values <- prediction_values %>% group_by(gwid,year) %>% arrange(gwid,year) %>% mutate(tek_state = median(tek_state,na.rm=T))#PS institutions targeting each group within a given country year on average
prediction_values <- unique(prediction_values[c("country","gwid","year", "geo_conc", "d10_victory_neg", "tek_state", "minoritysum","ongoing_grp","loggdppc","logpop","democracy")])

##predicting the effect of corporate power-sharing
Belgium <- Predict(m5.1, ps1h_corp_g_strength = seq(0,1,by=0.0050251), mminority = c(0,1), ps1h_corp_mo_others = 0, ps1h_lib_mo = 0, geo_conc = 1, tek_state = 1, minoritysum = 0.4100000, democracy = 1, gwid = c(211), year = 2013, ongoing_grp = 0, d10_victory_neg = 0, loggdppc = 10.619573, logpop = 2.4143784, conf.int = 0.9, adj.zero = T, fun=plogis, kint = 2)
Switzerland <- Predict(m5.1, ps1h_corp_g_strength = seq(0,1,by=0.0050251), mminority = c(0,1), ps1h_corp_mo_others = 0, ps1h_lib_mo = 0, geo_conc = 1, tek_state = 1, democracy = 1, gwid = c(225), year = 2013, ongoing_grp = 0, d10_victory_neg = 0, minoritysum = 0.3780000, loggdppc = 10.942628, logpop = 2.0905479, conf.int = 0.9, adj.zero = T, fun=plogis, kint = 2)
South_Africa <- Predict(m5.1, ps1h_corp_g_strength = seq(0,1,by=0.0050251), mminority = c(0,1), ps1h_corp_mo_others = 0, ps1h_lib_mo = 0, geo_conc = 1, tek_state = 0, democracy = 0, gwid = c(560), year = 2013, ongoing_grp = 0, d10_victory_neg = 0, minoritysum = 0.7700000, loggdppc = 9.429125, logpop = 3.9761606, conf.int = 0.9, adj.zero = T, fun=plogis, kint = 2)
Bosnia <- Predict(m5.1, ps1h_corp_g_strength = seq(0,1,by=0.0050251), mminority = c(0,1), ps1h_corp_mo_others = 0, ps1h_lib_mo = 0, geo_conc = 1, tek_state = 1, democracy = 0, gwid = c(346), year = 2013, ongoing_grp = 0, d10_victory_neg = 0, minoritysum = 0.5160000, loggdppc = 9.238405, logpop = 1.2823215, conf.int = 0.9, adj.zero = T, fun=plogis, kint = 2)
Kosovo <- Predict(m5.1, ps1h_corp_g_strength = seq(0,1,by=0.0050251), mminority = c(0,1), ps1h_corp_mo_others = 0, ps1h_lib_mo = 0, geo_conc = 0, tek_state = 1, democracy = 1, gwid = c(347), year = 2013, ongoing_grp = 0, d10_victory_neg = 0, minoritysum = 0.1000000, loggdppc = 9.051347, logpop = 0.6010867, conf.int = 0.9, adj.zero = T, fun=plogis, kint = 2)
Macedonia <- Predict(m5.1, ps1h_corp_g_strength = seq(0,1,by=0.0050251), mminority = c(0,1), ps1h_corp_mo_others = 0, ps1h_lib_mo = 0, geo_conc = 0, tek_state = 1, democracy = 1, gwid = c(343), year = 2013, ongoing_grp = 0, d10_victory_neg = 0, minoritysum = 0.3580000, loggdppc = 9.382415, logpop = 0.7303172, conf.int = 0.9, adj.zero = T, fun=plogis, kint = 2)
Nepal <- Predict(m5.1, ps1h_corp_g_strength = seq(0,1,by=0.0050251), mminority = c(0,1), ps1h_corp_mo_others = 0, ps1h_lib_mo = 0, geo_conc = 1, tek_state = 0, democracy = 1, gwid = c(790), year = 2013, ongoing_grp = 0, d10_victory_neg = 1, minoritysum = 0.3800000, loggdppc = 7.679691, logpop = 3.3316797, conf.int = 0.9, adj.zero = T, fun=plogis, kint = 2)
Lebanon <- Predict(m5.1, ps1h_corp_g_strength = seq(0,1,by=0.0050251), mminority = c(0,1), ps1h_corp_mo_others = 0, ps1h_lib_mo = 0, geo_conc = 1, tek_state = 0, democracy = 0, gwid = c(660), year = 2013, ongoing_grp = 0, d10_victory_neg = 0, minoritysum = 0.6800000, loggdppc = 9.575028, logpop = 1.6631876, conf.int = 0.9, adj.zero = T, fun=plogis, kint = 2)
Nigeria <- Predict(m5.1, ps1h_corp_g_strength = seq(0,1,by=0.0050251), mminority = c(0,1), ps1h_corp_mo_others = 0, ps1h_lib_mo = 0, geo_conc = 1, tek_state = 0, democracy = 1, gwid = c(475), year = 2013, ongoing_grp = 0, d10_victory_neg = 0, minoritysum = 0.7100000, loggdppc = 8.608689, logpop = 5.1465016, conf.int = 0.9, adj.zero = T, fun=plogis, kint = 2)
Malaysia <- Predict(m5.1, ps1h_corp_g_strength = seq(0,1,by=0.0050251), mminority = c(0,1), ps1h_corp_mo_others = 0, ps1h_lib_mo = 0, geo_conc = 1, tek_state = 0, democracy = 0, gwid = c(820), year = 2013, ongoing_grp = 0, d10_victory_neg = 0, minoritysum = 0.4230000, loggdppc = 10.052950, logpop = 3.3913734, conf.int = 0.9, adj.zero = T, fun=plogis, kint = 2)
predictions <- rbind("Belgium" = Belgium, "Switzerland" = Switzerland, "Bosnia" = Bosnia, "Kosovo" = Kosovo, "Macedonia" = Macedonia, "Lebanon" = Lebanon, "Malaysia" = Malaysia, "Nepal" = Nepal, "South Africa" = South_Africa, "Nigeria" = Nigeria)

##predicting the effect of corporate power-sharing
Belgium2 <- Predict(m5.1, ps1h_corp_g_strength = seq(0,1,by=0.0050251), mminority = c(0,1), ps1h_corp_mo_others = 1, ps1h_lib_mo = 0, geo_conc = 1, tek_state = 1, minoritysum = 0.4100000, democracy = 1, gwid = c(211), year = 2013, ongoing_grp = 0, d10_victory_neg = 0, loggdppc = 10.619573, logpop = 2.4143784, conf.int = 0.9, adj.zero = T, fun=plogis, kint = 2)
Switzerland2 <- Predict(m5.1, ps1h_corp_g_strength = seq(0,1,by=0.0050251), mminority = c(0,1), ps1h_corp_mo_others = 1, ps1h_lib_mo = 0, geo_conc = 1, tek_state = 1, democracy = 1, gwid = c(225), year = 2013, ongoing_grp = 0, d10_victory_neg = 0, minoritysum = 0.3780000, loggdppc = 10.942628, logpop = 2.0905479, conf.int = 0.9, adj.zero = T, fun=plogis, kint = 2)
South_Africa2 <- Predict(m5.1, ps1h_corp_g_strength = seq(0,1,by=0.0050251), mminority = c(0,1), ps1h_corp_mo_others = 1, ps1h_lib_mo = 0, geo_conc = 1, tek_state = 0, democracy = 0, gwid = c(560), year = 2013, ongoing_grp = 0, d10_victory_neg = 0, minoritysum = 0.7700000, loggdppc = 9.429125, logpop = 3.9761606, conf.int = 0.9, adj.zero = T, fun=plogis, kint = 2)
Bosnia2 <- Predict(m5.1, ps1h_corp_g_strength = seq(0,1,by=0.0050251), mminority = c(0,1), ps1h_corp_mo_others = 1, ps1h_lib_mo = 0, geo_conc = 1, tek_state = 1, democracy = 0, gwid = c(346), year = 2013, ongoing_grp = 0, d10_victory_neg = 0, minoritysum = 0.5160000, loggdppc = 9.238405, logpop = 1.2823215, conf.int = 0.9, adj.zero = T, fun=plogis, kint = 2)
Kosovo2 <- Predict(m5.1, ps1h_corp_g_strength = seq(0,1,by=0.0050251), mminority = c(0,1), ps1h_corp_mo_others = 1, ps1h_lib_mo = 0, geo_conc = 0, tek_state = 1, democracy = 1, gwid = c(347), year = 2013, ongoing_grp = 0, d10_victory_neg = 0, minoritysum = 0.1000000, loggdppc = 9.051347, logpop = 0.6010867, conf.int = 0.9, adj.zero = T, fun=plogis, kint = 2)
Macedonia2 <- Predict(m5.1, ps1h_corp_g_strength = seq(0,1,by=0.0050251), mminority = c(0,1), ps1h_corp_mo_others = 1, ps1h_lib_mo = 0, geo_conc = 0, tek_state = 1, democracy = 1, gwid = c(343), year = 2013, ongoing_grp = 0, d10_victory_neg = 0, minoritysum = 0.3580000, loggdppc = 9.382415, logpop = 0.7303172, conf.int = 0.9, adj.zero = T, fun=plogis, kint = 2)
Nepal2 <- Predict(m5.1, ps1h_corp_g_strength = seq(0,1,by=0.0050251), mminority = c(0,1), ps1h_corp_mo_others = 1, ps1h_lib_mo = 0, geo_conc = 1, tek_state = 0, democracy = 1, gwid = c(790), year = 2013, ongoing_grp = 0, d10_victory_neg = 1, minoritysum = 0.3800000, loggdppc = 7.679691, logpop = 3.3316797, conf.int = 0.9, adj.zero = T, fun=plogis, kint = 2)
Lebanon2 <- Predict(m5.1, ps1h_corp_g_strength = seq(0,1,by=0.0050251), mminority = c(0,1), ps1h_corp_mo_others = 1, ps1h_lib_mo = 0, geo_conc = 1, tek_state = 0, democracy = 0, gwid = c(660), year = 2013, ongoing_grp = 0, d10_victory_neg = 0, minoritysum = 0.6800000, loggdppc = 9.575028, logpop = 1.6631876, conf.int = 0.9, adj.zero = T, fun=plogis, kint = 2)
Nigeria2 <- Predict(m5.1, ps1h_corp_g_strength = seq(0,1,by=0.0050251), mminority = c(0,1), ps1h_corp_mo_others = 1, ps1h_lib_mo = 0, geo_conc = 1, tek_state = 0, democracy = 1, gwid = c(475), year = 2013, ongoing_grp = 0, d10_victory_neg = 0, minoritysum = 0.7100000, loggdppc = 8.608689, logpop = 5.1465016, conf.int = 0.9, adj.zero = T, fun=plogis, kint = 2)
Malaysia2 <- Predict(m5.1, ps1h_corp_g_strength = seq(0,1,by=0.0050251), mminority = c(0,1), ps1h_corp_mo_others = 1, ps1h_lib_mo = 0, geo_conc = 1, tek_state = 0, democracy = 0, gwid = c(820), year = 2013, ongoing_grp = 0, d10_victory_neg = 0, minoritysum = 0.4230000, loggdppc = 10.052950, logpop = 3.3913734, conf.int = 0.9, adj.zero = T, fun=plogis, kint = 2)
predictions2 <- rbind("Belgium" = Belgium2, "Switzerland" = Switzerland2, "Bosnia" = Bosnia2, "Kosovo" = Kosovo2, "Macedonia" = Macedonia2, "Lebanon" = Lebanon2, "Malaysia" = Malaysia2, "Nepal" = Nepal2, "South Africa" = South_Africa2, "Nigeria" = Nigeria2)

##predicting the effect of liberal power-sharing##
Belgium3 <- Predict(m5.1, ps1h_lib_mo = seq(0,1,by=0.0050251), mminority = c(0,1), ps1h_corp_mo_others = 0, ps1h_corp_g_strength = 0, geo_conc = 1, tek_state = 1, democracy = 1, gwid = c(211), year = 2013, ongoing_grp = 0, d10_victory_neg = 0, minoritysum = 0.4100000, loggdppc = 10.619573, logpop = 2.4143784, conf.int = 0.9, adj.zero = T, fun=plogis, kint = 2)
Switzerland3 <- Predict(m5.1, ps1h_lib_mo = seq(0,1,by=0.0050251), mminority = c(0,1), ps1h_corp_mo_others = 0, ps1h_corp_g_strength = 0, geo_conc = 1, tek_state = 1, democracy = 1, gwid = c(225), year = 2013, ongoing_grp = 0, d10_victory_neg = 0, minoritysum = 0.3780000, loggdppc = 10.942628, logpop = 2.0905479, conf.int = 0.9, adj.zero = T, fun=plogis, kint = 2)
South_Africa3 <- Predict(m5.1, ps1h_lib_mo = seq(0,1,by=0.0050251), mminority = c(0,1), ps1h_corp_mo_others = 0, ps1h_corp_g_strength = 0, geo_conc = 1, tek_state = 0, democracy = 0, gwid = c(560), year = 2013, ongoing_grp = 0, d10_victory_neg = 0, minoritysum = 0.7700000, loggdppc = 9.429125, logpop = 3.9761606, conf.int = 0.9, adj.zero = T, fun=plogis, kint = 2)
Bosnia3 <- Predict(m5.1, ps1h_lib_mo = seq(0,1,by=0.0050251), mminority = c(0,1), ps1h_corp_mo_others = 0, ps1h_corp_g_strength = 0, geo_conc = 1, tek_state = 1, democracy = 0, gwid = c(346), year = 2013, ongoing_grp = 0, d10_victory_neg = 0, minoritysum = 0.5160000, loggdppc = 9.238405, logpop = 1.2823215, conf.int = 0.9, adj.zero = T, fun=plogis, kint = 2)
Kosovo3 <- Predict(m5.1, ps1h_lib_mo = seq(0,1,by=0.0050251), mminority = c(0,1), ps1h_corp_mo_others = 0, ps1h_corp_g_strength = 0, geo_conc = 0, tek_state = 1, democracy = 1, gwid = c(347), year = 2013, ongoing_grp = 0, d10_victory_neg = 0, minoritysum = 0.1000000, loggdppc = 9.051347, logpop = 0.6010867, conf.int = 0.9, adj.zero = T, fun=plogis, kint = 2)
Macedonia3 <- Predict(m5.1, ps1h_lib_mo = seq(0,1,by=0.0050251), mminority = c(0,1), ps1h_corp_mo_others = 0, ps1h_corp_g_strength = 0, geo_conc = 0, tek_state = 1, democracy = 1, gwid = c(343), year = 2013, ongoing_grp = 0, d10_victory_neg = 0, minoritysum = 0.3580000, loggdppc = 9.382415, logpop = 0.7303172, conf.int = 0.9, adj.zero = T, fun=plogis, kint = 2)
Nepal3 <- Predict(m5.1, ps1h_lib_mo = seq(0,1,by=0.0050251), mminority = c(0,1), ps1h_corp_mo_others = 0, ps1h_corp_g_strength = 0, geo_conc = 1, tek_state = 0, democracy = 1, gwid = c(790), year = 2013, ongoing_grp = 0, d10_victory_neg = 1, minoritysum = 0.3800000, loggdppc = 7.679691, logpop = 3.3316797, conf.int = 0.9, adj.zero = T, fun=plogis, kint = 2)
Lebanon3 <- Predict(m5.1, ps1h_lib_mo = seq(0,1,by=0.0050251), mminority = c(0,1), ps1h_corp_mo_others = 0, ps1h_corp_g_strength = 0, geo_conc = 1, tek_state = 0, democracy = 0, gwid = c(660), year = 2013, ongoing_grp = 0, d10_victory_neg = 0, minoritysum = 0.6800000, loggdppc = 9.575028, logpop = 1.6631876, conf.int = 0.9, adj.zero = T, fun=plogis, kint = 2)
Nigeria3 <- Predict(m5.1, ps1h_lib_mo = seq(0,1,by=0.0050251), mminority = c(0,1), ps1h_corp_mo_others = 0, ps1h_corp_g_strength = 0, geo_conc = 1, tek_state = 0, democracy = 1, gwid = c(475), year = 2013, ongoing_grp = 0, d10_victory_neg = 0, minoritysum = 0.7100000, loggdppc = 8.608689, logpop = 5.1465016, conf.int = 0.9, adj.zero = T, fun=plogis, kint = 2)
Malaysia3 <- Predict(m5.1, ps1h_lib_mo = seq(0,1,by=0.0050251), mminority = c(0,1), ps1h_corp_mo_others = 0, ps1h_corp_g_strength = 0, geo_conc = 1, tek_state = 0, democracy = 0, gwid = c(820), year = 2013, ongoing_grp = 0, d10_victory_neg = 0, minoritysum = 0.4230000, loggdppc = 10.052950, logpop = 3.3913734, conf.int = 0.9, adj.zero = T, fun=plogis, kint = 2)
predictions3 <- rbind("Belgium" = Belgium3, "Switzerland" = Switzerland3, "Bosnia" = Bosnia3, "Kosovo" = Kosovo3, "Macedonia" = Macedonia3, "Lebanon" = Lebanon3, "Malaysia" = Malaysia3, "Nepal" = Nepal3, "South Africa" = South_Africa3, "Nigeria" = Nigeria3)

##figure 4a: predicted probabilities, corporate power-sharing (no corporate provisions for other groups)
figure4a <- ggplot(predictions, ylab = "De facto political status", xlab = "Strength of de jure corporate power-sharing (group)", legend.label = "Micro-minority") + ggtitle("a) Corporate power-sharing (no corporate provisions targeting other groups)") + scale_x_continuous(limits = c(0, 1)) +
  theme(legend.position = "bottom", axis.text = element_text(size=5), legend.text = element_text(size = 6), plot.title = element_text(size=10), text=element_text(family="Times New Roman", size = 8))

##figure 4b: predicted probabilities, corporate power-sharing (full corporate provisions for other groups)
figure4b <- ggplot(predictions2, ylab = "De facto political status", xlab = "Strength of de jure corporate power-sharing (group)", legend.label = "Micro-minority") + ggtitle("b) Corporate power-sharing (full corporate provisions targeting other groups)") + scale_x_continuous(limits = c(0, 1))  +
  theme(legend.position = "bottom", axis.text = element_text(size=5), legend.text = element_text(size = 6), plot.title = element_text(size=10), text=element_text(family="Times New Roman", size = 8))

##figure 4c: predicted probabilities, liberal power-sharing (no corporate provisions for other groups)
figure4c <- ggplot(predictions3, ylab = "De facto political status", xlab = "Strength of de jure liberal power-sharing (country)", legend.label = "Micro-minority") + ggtitle("c) Liberal power-sharing (no corporate provisions targeting other groups)") + scale_x_continuous(limits = c(0, 1))  +
  theme(legend.position = "bottom", axis.text = element_text(size=5), legend.text = element_text(size = 6), plot.title = element_text(size=10), text=element_text(family="Times New Roman", size = 8))

##arranging them in one figure
library(gridExtra)
library(grid)
grid_arrange_shared_legend <- function(..., ncol = length(list(...)), nrow = 1, position = c("bottom", "right")) {
  plots <- list(...)
  position <- match.arg(position)
  g <- ggplotGrob(plots[[1]] + 
                    theme(legend.position = position))$grobs
  legend <- g[[which(sapply(g, function(x) x$name) == "guide-box")]]
  lheight <- sum(legend$height)
  lwidth <- sum(legend$width)
  gl <- lapply(plots, function(x) x +
                 theme(legend.position = "none"))
  gl <- c(gl, ncol = ncol, nrow = nrow)
  
  combined <- switch(position,
                     "bottom" = arrangeGrob(do.call(arrangeGrob, gl), 
                                            legend,ncol = 1,
                                            heights = unit.c(unit(1, "npc") - lheight, lheight)),
                     "right" = arrangeGrob(do.call(arrangeGrob, gl),
                                           legend, ncol = 2,
                                           widths = unit.c(unit(1, "npc") - lwidth, lwidth)))
  
  grid.newpage()
  grid.draw(combined)
  
  # return gtable invisibly
  invisible(combined)
}#function to arrange plots
figure4 <- grid_arrange_shared_legend(figure4a,figure4b, figure4c, ncol=1,nrow=3)
ggsave(file="figure4.jpg", figure4, width = 20.1, height = 13.8, units ="cm", dpi = 800)
